Characterization of the gravitational wave emission of three black holes 



Pablo Galaviz 1, 2 [] and Bernd Briigmann 2 '[j] 

1 School of Mathematical Science, Monash University, Melbourne, VIC 3800, Australia 
2 Theoretical Physics Institute, University of Jena, 07743 Jena, Germany 
(Dated: August 24, 2011) 

We study the gravitational wave emission of three compact objects using post-Newtonian (PN) 
equations of motion derived from the Arnowitt-Deser-Misner Hamiltonian formulation, where we in- 
clude (for the first time in this context) terms up to 2.5 PN order. We perform numerical simulations 
of a hierarchical configuration of three compact bodies in which a binary system is perturbed by a 
third, lighter body initially far away from the binary. The relative importance of the different PN 
orders is examined. We compute the waveform in the linear regime considering mass quadrupole, 
current quadrupole and mass octupole contributions. Performing a spherical harmonic decomposi- 
tion of the waveforms we find that from the / = 3 modes it is possible to extract information about 
the third body, in particular, the period, eccentricity of its orbit, and the inclination angle between 
the inner and outer binary orbits. 
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I. INTRODUCTION 

In the near future gravitational wave detectors will 
open a new branch on astronomy beyond electromag- 
netism and particles. In order to extract astrophysical 
information from the waves it is crucial to model the 
source in an accurate way. One of the parameters to 
determine is the number of bodies which generate the 
waves. Some configurations of three bodies can produce 
particular periodic waveforms with distinctive features, 
e.g. pQ. On the other hand, Lagrange's triangle solu- 
tion produces a quadrupole waveform which is identical 
to the one produced by a binary system [2J. However, it 
is possible to distinguish between a binary system and a 
triple one by considering the octupole part of the wave- 
form [3J. Lagrange's solution is stable only if one of the 
bodies holds more than 95% of the total mass [4]. For 
three-body systems of comparable masses there are sta- 
ble configuration in which it is possible to characterize 
the system by looking at the waveform. In this work 
we consider Jacobian systems, in which the three-body 
configuration is composed of two parts, a clearly defined 
binary and a third body orbiting far away. We will refer 
to this kind of system also as a hierarchical system. 

Several models of three or more black holes were re- 
cently studied from the astrophysical point of view. Hi- 
erarchical three black hole configurations interacting in 
a galactic core were studied by several authors. For ex- 
ample in [MI] some configurations of intermediate-mass 
black holes with different mass ratios were considered. 
The inclusion of gravitational radiation was done via an 
effective force which includes 1 PN (post-Newtonian) and 
2.5 PN corrections to the binary dynamics. The config- 
urations consist of a binary system in a quasicircular or- 
bit and a third black hole approaching from a distance 
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around 200 times the binary separation. The initial ec- 
centricity was specified in a random way. N-body simu- 
lations of dynamical evolution of triple equal-mass super- 
massive black holes in a galactic nuclei were performed 
in [8] . The method includes an effective force with gravi- 
tational radiation terms and galaxy halo interactions. In 
[5] the dynamics of repeated triple supermassive black 
hole interactions in galactic nuclei with several mass ra- 
tios and eccentricities were considered. The simulations 
were performed using Newtonian dynamics with correc- 
tions through an additional force which includes 2.5 PN 
corrections to the binary dynamics and stellar dynami- 
cal friction. Other astrophysical applications of multiple 
black hole simulations include, for example, three-body 
kicks [ini E] and binary-binary encounters (see e.g. [T2T - 

USD- 

The first complete simulations using general- 
relativistic numerical evolutions of three black holes 
were presented in [T71 [18]. These recent simulations 
show that the dynamics of three compact objects display 
a qualitatively different behavior than the Newtonian 
dynamics. In [TH] the sensitivity of fully relativistic 
evolutions of three and four black holes to changes in 
the initial data was examined, where the examples for 
three black holes are some of the simpler cases already 
discussed in [TTIIH]- The apparent horizon and the event 
horizon of multiple black holes have been studied in 
[2"Ull22) . Although fully general-relativistic simulations 
are available, they are limited to only a small number of 
orbits for small separations of the black holes. 

In the present work we study three-body systems with 
PN methods, where the main technical novelty is the in- 
clusion of the 2.5 PN terms in the orbital dynamics. We 
do not consider compact objects with spin, although re- 
cently the knowledge of Hamiltonian up to 2.5 PN was 
completed with the computation of a next-to-leading or- 
der spin-orbit and spin-spin Hamiltonian 23 25 

Using post-Newtonian techniques, it is currently possi- 
ble to describe the dynamics of n compact objects with- 
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out spin up to 3.5 PN order (see e.g. [23H32]), although 
explicit and closed expressions for the terms required for 
our purpose are not available for arbitrary n. For binary 
systems the Arnowitt-Deser-Misner (ADM) Hamiltonian 
has been specialized up to 3.5 PN order [S3]. For three 
bodies there are explicit formulas up to 2.5 PN order, 
with the final key integral performed in [33] (see also 
Appendix [A] and [35] ; see [35] for an alternative discus- 
sion where the final result is not as yet explicit) . For four 
compact objects, the same degree of explicitness has not 
been obtained, see e.g. [37] on which [33] is based where 
the integral (C.2) must be computed for four bodies, and 
also Appendix D of [35]. In the ADM-Hamiltonian for- 
malism the resulting equations of motion exactly con- 
serve the constants of motion. For numerical simulations 
this represents an advantage with respect to other post- 
Newtonian approaches, since the constants of motion can 
be tracked and their constancy continually checked. 

Periodic solutions were studied using the 1 PN and 
2 PN approximation in (35j [39j 00] . Examples of three 
compact bodies in a collinear configuration were consid- 
ered in [3UI32], and Lagrange's equilateral triangular so- 
lution was studied including 1 PN effects in [33] . In [33] , 
the stability of the Lagrangian points in a black hole bi- 
nary system was studied in the test particle limit, for 
which radiation effects were modeled by a drag force. 

The most likely source of gravitational waves are bi- 
nary compact objects. Recently it was shown that the 
probability that more than two black holes interact in 
the strongly relativistic regime is, not surprisingly, very 
small [45] . For practical purposes the creation of gravita- 
tional waveform templates for gravitational wave detec- 
tors is naturally focused on binary systems, and even bi- 
nary systems can produce complicated waveforms when 
taking into account spinning black holes and eccentric 
orbits, e.g. @B]. 

Nevertheless, it remains an interesting question of prin- 
ciple what additional wave phenomena are possible for 
more than two compact objects. The waveform char- 
acterization of three or more compact objects is comple- 
mentary to the study of binary systems. For a three-body 
system the complexity of the orbits can reveal properties 
of the waves which for a binary system are hidden. As 
we will demonstrate for a hierarchical system, from the 
I = 3 modes of the gravitational wave it is possible to ex- 
tract information about the third body, particularly the 
period, eccentricity of its orbit and the inclination angle 
between the inner and outer binary orbits. 

The paper is organized as follows. In Sec. [TTJ we sum- 
marize the equation of motion up to 2.5 post-Newtonian 
approximation for three bodies. This is followed by a 
discussion of gravitational radiation in the linear regime 
and the multipole expansion of the gravitational waves. 
In Sec. |III A| we describe the numerical techniques used 
to solve the equation of motion and we present some re- 
sults for test cases. The perturbation of a binary system 
by a third object is presented in Sec. |III B] where we per- 
form numerical experiments in order to characterize the 



waveform. We conclude in Sec. HV 



A. Notation and units 

We employ the following notation: x = (xi) denotes a 
point in the three-dimensional Euclidean space M 3 , let- 
ters a,b, . . . are particle labels. We define r a := x — x a , 
r a := |r a |, n a := r a /r a ; for a ^ b, r ab := x a - x b , 
Tab ■= \f a b\, n ab := f ab /r ab ; here | • | denotes the length 
of a vector. The mass parameter of the a-th particle 
is denoted by m a , with M = X)a m a- Summation runs 
from 1 to 3. The linear momentum vector is denoted by 
p a - A dot over a symbol, as in x, means the total time 
derivative, and partial differentiation with respect to x l 
is denoted by 0j. 

In order to simplify the calculations it is useful to define 
dimensionless variables (see e.g. [27]). We use as basis 
quantities for the Newtonian and post-Newtonian calcu- 
lation the gravitational constant G, the speed of light c 
and the total mass of the system M. Using derived con- 
stants for time r = A/G/c 3 , length I — MG/c 2 : linear 
momentum V = Mc and energy E = Mc 2 we construct 
dimensionless variables. The physical variables are re- 
lated with the dimensionless variables by means of a scal- 
ing, for example, denoting with capital letters the phys- 
ical variables with the usual dimensions and with lower- 
case the dimensionless variable we define for a particle a 
its position x a := X a /l, linear momentum p a := P a /V 
and mass m a = M a /M (notice that m a < 1, Va). 



II. EQUATIONS OF MOTION 

In the ADM post-Newtonian approach it is possible to 
split the Hamiltonian in a series with coefficients which 
are inverse powers of the speed of light (see e.g. [2T1I25] ) 



u 



<2.5 



l n 2 + c~ 5 h 



2.5- 



(1) 



Here each term of the Hamiltonian c n H n /2 is a quantity 
with a dimension of energy, and we write it explicitly 
with factors of c. The dimensionless Hamiltonian is given 
by H n / 2 — c n 'H n /2/S. For each term we calculate the 
equations of motion 



(4)n 
-(Pa)n 



dHn 

dpi ' 

dHn 

dxi ' 



(2) 
(3) 



where the equations of motion up to 2.5 PN approxima- 
tion are 

X a = (x a )a + (x a )i + (Xa)2 + (x a )2.5, ( 4 ) 
Pa = (Pa)o + (Pa)l + (Pa)2 + {Pah.h- (5) 
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The first term in ([!]) is the Hamiltonian for n particles 
interacting under Newtonian gravity, 



H 



1 n 2 

z \ " Pa 

2 ^ m„ 



n 

a.b^a 



m a m b 

Tab 



(6) 



with p a — m a x a 2 . The inclusion of post-Newtonian cor- 
rections enriches the phenomenology of the system. 



1. Post- Newtonian equations of motion up to 2.5 order 

The first post-Newtonian correction to the equations of 
motion is discussed extensively in the literature (see e.g. 
[271 0H] ) • The three-body Hamiltonian at first and sec- 
ond post- Newtonian order is given in Appendix [A] The 
equations of motion for the first post-Newtonian order 
are given by pi (p§ and (pUL For particle a we obtain 



Pa 

2m? 



'« - \ ^ ~~~ \6—Pa - 7Pb - {flab ' Pbjflab 



\ m a 

bjta 



(7) 



(Pa)l = -9^ 



b^a 



3— Pa + 3—Pb ~ 7 (Pa ■ Pb) - 3{fl a b ' Pa){hab ' Pa) 

m a m b 



^2 



EE 



m a m b m c 



hab + E E 



m a m b m c 



1 



2 '"ao 1 / / 2 Tlab 

r ab r ac h=Ln r=Lh r ab r bc 



E 

a^b 



(«ab ' Pb)Pa + (n a b ' Pa)Pb 



(8) 



For the second post-Newtonian approximation the equations of motion are calculated using (|]), ([3| and (|A2]). For 
brevity we do not display the explicit equations. 

Following [33] we obtain equations of motion from the 2.5 PN Hamiltonian in the ADM gauge. The general 
2.5 PN Hamiltonian is 



#2.5 = J^X(4)ij{Xa',Pa']t)X(4)ij{Xa,Pa), 



(9) 



where the auxiliary function X(4)ij i s defined by 



X(4)ij(Xa,Pa) ■■= — (Pa 6 ij ~ ^PaiPaj) + E E 



m a m b 

Tab 



(10) 



a a b ^La 

Our expressions differ from [261 133] due to a different choice of units. The explicit form of the derivative in (J9J) is 

2 

m„ 



2 r 

X(4)ij(x a >,Pa>) = ^>^ \2(Pa> ' Pa')5ij ~ 3(pa>iPa>j + Pa'iPa'j) 

if in 1 L 



+EE 



m a >m b , 



a' b'=jta' Ta ' b ' 



[3(r a 'b'in a 'b'j + na'b'ifa'b'j) + {n a >v •fVb'X^j - 9n a >b'in a 'b'j) 



(11) 



We denote the retarded variables by primed quantities. The position and momentum appearing in Eq. (Ill are 
not affected by the derivative operators given by ([2| and ([3]), and only after calculating those derivatives we identify 
positions and momenta inside and outside the transverse-traceless variables (i.e. the primed and unprimed quantities). 
We replace the time derivatives of the primed coordinates and positions given in Eq. (11) by the 1 PN equations of 
motion Eqs. (7]) and 

The equations of motion for 2.5 PN are given in short hand by 



1 d 

(X a )2.5 = ^X(4)ij{Xa,Pa](x a ) 1 ,(p a )l,t)—X(4)ij(Xa,Pa), 

1 d 

(Pa)2.5 = --^X(4,)ij{X a ,p a ; (X a )i,(p )i, t)—X(i)ij(^a, Pa)- 



(12) 
(13) 



Given initial values for x a and p a of each particle it numerically, 
is possible to integrate the resulting equations of motion 
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A. Gravitational radiation in the linear regime 

We consider leading order and next-to-leading order 
gravitational waves calculated using trajectories which 
contain post-Newtonian corrections. We compute the 
gravitational waveforms for a given observational direc- 
tion, and alternatively we calculate the multipole decom- 
position which allows us to reconstruct the waves for an 
arbitrary direction. The inclusion of post-Newtonian cor- 
rections to the gravitational waveforms is a topic for fu- 
ture research in the thrce-compact-body problem. 



1. Quadrupole and octupole formulas 

Here we summarize the formulas for quadrupole and 
octupole mass radiation and for current quadrupole ra- 
diation (for a review see e.g. [HI [SO])- The second and 
third mass moments arc defined by 



M«(i)= / T Q0 (x,t)x l x ] d 3 x, 
The second moment of the momentum density is 
For n point particles 



(14) 
(15) 

(16) 
(17) 



where 7 a := (1 — p a 2 ) 1 / 2 is the Lorentz factor, and p% 



Ja{ m aiPa) is the four-momentum. In this case Eqs (14) 
( 16 1 reduce to 



M^(t)=J2^m a xi(t)xi(t), (18) 

a 

M llk {t) =Y,^m a x l a (t)xi(t)x k a (t), (19) 

a 

P^ k {t)=Y,^t)<{t)<{t)- (20) 

a 

In the following we consider the case where \p a \ <C 1, 
The mass quadrupole and octupole moment are given 

by 

Q^it) = M lj ~ \s i3 M kk , 

O ljk {t) = M tjk - -(5 ij M Uk + 5 lk M l]l + S jk M M ), 
5 

where repeated indices mean summation from 1 to 3. The 
current quadrupole is given by 



C k ' lm (t) = P 



k.lm 



p 



I, krn 



2P 



rn.kl 



A projection tensor into the plane nor- 
mal to the direction of wave propagation, 
h = (sin 8 sin <j>, sin 9 cos </>, cos 8), is defined by 



Ayfcz(n) := V ik Vji - -VijVki- 



(22) 
(23) 



The mass quadrupole and octupole waveforms are given 
by 



hJj r (x,t) MQ = r A ijkl (n)Q kl (t-r), 



(24) 



fifi r (x, t) MO = ^^ 3 ki^) n mO klm (t - r), (25) 



3r 



and the current quadrupole contribution to the waveform 
is 



k.ll: 



h ij {x,t) C Q = 7^^ijki(n)n m C ' 
The total contribution on the waveform is given by 

h[j T (x,t) =hJ J T (x,t) MQ + hJ ] T {x,t) CQ 
+ hf^(x,t) MO + .... 



(26) 



(27) 



where . . . means additional multipolcs. Assuming that 
the wave propagates in the z-direction, then h + = hj^f 
and h x = hfj ■ For an arbitrary direction n(8, <fi) we have 
to perform a rotation of the axes in order to identify the 
polarization with the hj^ and /i^ 2 T components. 

We decompose /i+ and h x into modes using spherical 
harmonics with spin- weight minus two, 



h + -zh x =J2 E -2^(e,$)hr, 



(28) 



/ m— — l 



where 



,^(0,$) := (-l) S /^ T ^(-)( e ) e4m$ ( 29 ) 



4,(0) : = f] { ~ m{l + m)l( - 1 ~ m) - {l + S)l{l ~ S)!]1/2 



t=Ci 



(l + m- t)\(l - s - t)W.(t + s-m)\ 



x (cos 9/2) 



2l+m-s-2t 



(sine/2) 2 



(30) 



with C\ — max(0,?7i — s) and C2 = min(^ + m,l — s). 
Using the orthonormality of the spherical harmonics it is 
possible to compute h. l m by 

h l m = r - a Yl(B,n/2-<l>)(h + -ih x )dn, (31) 
Jo Jo 



(21) where dfi = sin SdGAcp. 
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III. SIMULATIONS AND RESULTS 
A. Numerical integration 

We solved the equations of motion numerically using 
Mathematica 7.0 [51 J. We used the built-in low-level 
functions of the NDSolve routine with a "double-step" 
method using as subalgorithm the "explicit midpoint" 
method. We divided long simulations into substeps in or- 
der to store the result from time to time and to avoid sat- 
urating the random-access-memory. With this approach 
we can produce accurate numerical solutions of the equa- 
tions of motion. For our purpose the performance of 
Mathematica solving the ordinary differential equation 
system is not an issue (see the performance and accuracy 
tests at the end of this section). 

An important issue in the numerical integration of a 
three-body system arises when two of the bodies come 
very close to each other. Adaptive step size methods can 
automatically maintain the necessary accuracy to prop- 
erly resolve the orbits in the close interaction, but issues 
of efficiency arise. For the Newtonian system a number 
of techniques have been developed that address problems 
with accuracy and efficiency, see e.g. [5^H55] and refer- 
ences therein. For our PN evolutions, efficiency was not 
an issue, and furthermore the equations of motion are 
not valid for arbitrarily small separation anyway. What 
is of relevance here is a convenient criterion of when to 
stop the evolution. We monitor the absolute value of 
each conservative part of the Hamiltonian relative to 
the sum of the absolute values, 



H?° := 100 



\H, 



Hi 



Iff 2 



(32) 



We stop the simulation when the contribution of the first 
post-Newtonian correction is larger than 10%. 

In the remainder of this section we report on several 
tests that allow us to estimate the numerical errors. We 
use the Lagrangian equilateral triangle solution to com- 
pare the numerical with an analytical solution. In La- 
grange's solution each body is sitting in one corner of 
an equilateral triangle (see e.g. [57]). We set the side of 
such triangle to L = 1000, the mass ratio to 1:2:3, and 
the eccentricity to zero. Then each body follows a circu- 
lar orbit (with different radii) around the center of mass. 
The solution in this case is not stable [4] , however for cir- 
cular orbits we can compute the waveforms and compare 
with the analytical expressions [3]. 

In Fig. [I] we show the relative variation of the Hamil- 
tonian 



AH := 



H(0) - H(t) 
H(0) ' 



(33) 



and for each body the relative variation of the position 
with respect to the center of mass. The variation of the 
Hamiltonian is small (close to machine accuracy), how- 
ever the error in the orbits grows fast, breaking the regu- 
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FIG. 1. Test using Lagrange's equilateral solution of the New- 
tonian three-body problem. Shown is the relative variation of 
the Hamiltonian (top) and the relative change in the orbits 
(bottom). 



[xlO-«] 



Error!/;'", 



8 
6 

i 4 

2 


8 
6 

x 

S3 4 

< 



—\ 1 1 - 



0.4 
0.2 




0.3 



0.6 



0.9 



_ i i i 



_l i i i L 



MO 
MO 
CQ 

(a) 



1 r 





1 1 1 1 1 1 1 1 1 1 1 1 1 


I 


. 0.4 






- 


- 0.2 

- 




L } , L( hi l> hi 

a' syJUy.' v"-."_,.v,v ^.\,v<vw WJif\ 


. — 


'- ( 


) 


. . i . . . i . . , i i , , i 1 

0.3 0.6 0.9 1.2 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



(b) 



0.3 



0.6 



0.9 



1.2 txltf] 



FIG. 2. Test using Lagrange's equilateral solution of the 
Newtonian three-body problem. Shown is the absolute value 
of the difference between the analytical expression and the 
numerical calculation for the mass quadrupole, mass octupole 
and the current quadrupole for each polarization of the wave- 
form. The insets show the mass octupole and the current 
quadrupole. 



lar trajectory. In this case, after seven orbits the numer- 
ical solution fails. The waves exhibit a similar behavior. 
In Fig. [2j we show the error for each polarization of the 
waveforms (l24|-(26). The error is defined as the absolute 



value of the difference between the numerical calculation 
and the analytical expression. The mass octupole ex- 
hibits a noisy error due to the complicated nature of the 
analytical expression. On the other hand, it seems that 
the error in the mass quadrupole starts growing before 
the errors in the mass octupole and current quadrupole. 
By looking at the analytical expressions this fact can be 
explained as follows (see Eqns. (B5|-(B10|). The mass 
quadrupole part contains a factor a 2 Lu z (where a is the 
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FIG. 3. Moore's figure eight solution. Relative variation of 
the Hamiltonian for a solution which includes 2 PN correc- 
tions. 



separation of the bodies and u> the orbital frequency). 
The mass octupole and current quadrupole have a fac- 



tor a 3 uj 3 . In terms of a the factors reduce to 



and 



a -3/2 respectively. A small change in the orbit is visi- 
ble at a smaller length scale, and then the growth in the 
waveforms seems to be delayed. 

We reproduce a few of the results from [35] , specifically 
the simulation of the equal-mass Moore's figure eight [55] . 
which includes first and second post-Newtonian correc- 
tions. Our choice of method was guided by numerical 
experiments to minimize the numerical error in this ex- 
ample. With the double-step, midpoint method we ob- 
tain fluctuations of the Hamiltonian of 10 -14 (see Fig. [3]), 
while other methods and parameter settings can show a 
significantly larger error. 

We tested our n— body 2.5 PN equations of motion for 
the case n — 2, i.e. for binary systems. The variation 
of the semimajor axis and of the eccentricity of a binary 
system due to the gravitational radiation is given by [58] 



da 

dT 

de 
dt 



64 



m 1 m 2 



5 a 3 (l-e 2 ) 7 /2 

304 77117712 

" 15 a 4 (l - e 2 ) 5 / 2 



1 H e 2 

24 

121 

f 304* 



37 L 
96 6 



(34) 
(35) 



We tested the 2.5 PN equations of motion (12) and (13) 



by comparison with direct numerical integration of the 
Eqs. (34) and (35). The test was performed with two 
different binaries, one with initial eccentricity eo = 0.1 
and one with eo = 0.5. In both cases we set m\ = 2m2, 
ao = 160. The numerical integration of the 2.5 PN equa- 
tions agree very well with the result provided by the nu- 
merical integration of (34) and (35). We calculate the 



eccentricity of our orbits with the Newtonian formula 



2l 2 H c 
(mim 2 ) 



(36) 



FIG. 4. Binary system with 2.5 PN radiation. Top: Relative 
variation of the apoapsis of the two bodies. Bottom: eccen- 
tricity variation; comparison of our numerical result (solid 
and dashed lines) with the numerical integration of ( 34 1 and 
(|35[) (marks + and x ) for two initial eccentricities. 
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FIG. 5. Henon criss-cross solution for Newtonian, 2 and 
2.5 PN dynamics. The main panel shows the relative variation 
of the Hamiltonian for the three cases. The inset shows only 
the conservative systems (Newtonian and 2 PN). 



where I is the magnitude of the total angular momen- 
tum and H c is the value of the conservative part of the 
Hamiltonian. The apoapsis (the maximum separation 
of the two bodies) is related to the semimajor axis by 
f*ap = a(l + e). For simplicity we compare in the upper 
panel of Fig. [4| the relative variation of r ap to its initial 
value and in the lower panel we show the variation of the 
eccentricity. 



In order to test the code for long evolutions of three 
bodies we use Henon's criss-cross solution [39l |59l ISO] . 
This solution is stable with respect to a wide range of 
perturbations [61] . We evolve the equal-mass criss-cross 
solution for around 10 3 orbits for ad-hoc initial parame- 
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FIG. 6. Henon criss-cross solution for Newtonian, 2 and 
2.5 PN dynamics. First and last orbits. In the Newtonian 
dynamics the orbits do not show a significant change. For 
dynamics including 2 PN corrections the orbits exhibit the 
expected precession. For the dynamics which includes 2.5 PN 
corrections the gravitational radiation produces a significant 
change in the orbits which in the long run breaks the system. 



FIG. 7. Hierarchical system. Initial configuration of the in- 
ner and external binaries. The initial momentum of the third 
body is given by considering the external binary as a Newto- 
nian binary. Shown are the osculating orbital planes IIi n and 
Ilext for inner and external binary orbits. The two planes are 
inclined by an angle i. 



ters. In our system of units, 

fi(0) = 1.07590A 2 i, pi(0) = 3~ 3 / 2 • O.^OdX^y, 
f 2 (0) = -0.07095A 2 x, p 2 (0) = -3~ 3 / 2 • 1.23187A- 1 y, 
x 3 (0) = -1.00496A 2 x, p 3 (0) = 3~ 3 / 2 • 1.03678A- 1 y, 

where x, y and z are the unitary basis vectors in Carte- 
sian coordinates, and A is a scaling factor (for our simula- 
tion A = 10). Notice that for this test we use the param- 
eters given in [61] with the scaling factor A, and doing a 
change of variables from initial velocity to initial momen- 
tum. Therefore, we are not including post-Newtonian 
corrections to the initial parameters. In Fig. [5] we plot 
the relative variation of the Hamiltonian for the evolu- 
tion using a Newtonian potential and the corresponding 
Hamiltonian variation for evolutions which include 2 and 
2.5 PN corrections. As is expected the variation of the 
Hamiltonian in the 2.5 PN case is huge compared to the 
conservative case, and the bodies separate after around 
t = 7.825 x 10 6 . The inner panel in Fig.[5]shows a detail 
of the conservative part. In this case the 2.5 PN dy- 
namics show better conservation of the Hamiltonian in 
contrast to the Newtonian case which has a variation in 
the Hamiltonian of around 4 x 10~ 12 . 

We confirm that the system is stable even after the 
inclusion of 2 and 2.5 PN corrections, see Fig. [6j In the 
Newtonian case the accumulation of numerical errors and 
probably a round-off in the initial parameters lead to a 
small variation of the orbits. The basic shape of the criss- 
cross figure suffers a small rotation. The 2 PN correction 
includes the effect of precession in the orbits; the original 
figure spins many times around the origin preserving its 
original shape. The inclusion of gravitational radiation 
via the 2.5 PN corrections has a stronger effect on the 
orbits, slowly deforming the original figure. The body in 



the circularlike orbit shows a significant reduction of the 
orbital radius, the two other bodies follow at the end a 
triangular orbit with narrow corners. 

We also use the Newtonian Henon criss-cross solution 
for performance and accuracy tests. A performance test 
based on walltime measurements resulted in about 4.4 
seconds per orbit on (one core of) an Intel i7-860 pro- 
cessor. For accuracy testing, we evaluate the error of 
time integration by a reversibility test. After computing 
a given number of orbits, we solve the system backward 
in time starting with the last position of each particle but 
replacing every linear momentum by its opposite value. 
To measure the error we compute the differences in phase 
space between the initial position and momentum and 
the position and momentum after the backward evolu- 
tion. For our standard setting, the error after 100 orbits 
is on the order of 10 -1 °. 



B. Strong perturbation of a binary system 

Here we consider the strong perturbation of the dy- 
namics and waveform of a binary compact object system 
due to a third smaller compact object. We take all PN 
corrections up to 2.5 PN for the three bodies. This ap- 
proach gives us a good description of the third body or- 
biting close to the binary. However, the computational 
cost of each simulation increases with respect to the New- 
tonian simulations, making it too costly to perform a 
comprehensive study of this study. Nevertheless, we can 
select a representative case in an attempt to identify key 
properties. 

As a basic configuration we study a Jacobian system 
with mass ratio 10:20:1. The inner binary system has 
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FIG. 8. Planar hierarchical system. Relative contribution to 
the Hamiltonian defined by ( 32 l . The inset shows the time 



dependence of H c — Ho + Hi + H2 . Notice that when the sys- 
tem approaches the merger phase, the Hamiltonian decreases 
quickly. 



initial separation r{,(0) = 150 and eccentricity e&(0) = 0. 
We set the initial parameters considering only the New- 
tonian dynamics, in particular, the eccentricity refers to 
the Newtonian case. We view the third compact body 
and the center of mass of the inner binary as a new bi- 
nary (we will refer to it as the external binary). The 
external binary has initial separation ^(0) = 10000 and 
initial eccentricity 63(0) = 0. The bodies start from a 
configuration where the apoapsis of the inner binary is 
perpendicular to the apoapsis of the external binary (see 
Fig.0. 

We denote the inclination angle between the osculat- 
ing orbital planes ITj n and n ext by i (see Fig. [7j. The 
behavior of the Hamiltonian is similar in every case that 
we consider. The conservative part of the Hamiltonian 
decreases relatively slowly during most of the simulation. 
However, when the system approaches the merger phase, 
the Hamiltonian decreases fast (see Fig. [8| . As we men- 
tioned before, the simulations are stopped when the con- 
tribution of the first post-Newtonian correction becomes 
larger than 10%. We consider this instant the time when 
the merger phase starts. 

We consider five numerical experiments. In Table [I] 
we summarize the configurations of the numerical exper- 
iments. We vary one parameter of the basic configuration 
and fix the rest. The main goal of the study is to charac- 
terize the changes produced in the waveforms due to the 
change in each parameter. 



1. Binary versus triple system 

We compare the case where the inner binary is not be- 
ing perturbed by the third compact body. Fig. [9] shows 
the components of the waveform for the h + polarization 
with an observational direction 9 = ir/4, = 0. In both 



TABLE I. Configuration of the numerical experiments. The 
fixed parameters in each case are the mass ratio 10 : 20 : 1, 
the initial eccentricity of the inner binary et = 0, and the an- 
gle between the apoapsis of the inner binary and the apoapsis 
of the external binary which is set to n/2. The base config- 
uration has initial binary separation r b = 150, Hamiltonian 
i/o+i+2+2.5, eccentricity of the external binary e 3 = 0, incli- 
nation angle of the osculating planes i — and initial external 
binary separation r 3 — 10000 



Experiment 


Parameter variation 


1 


r b G {130,140, 150,160,170} 


2 


H G {H +2.5, #0+1+2.5, -Ho+l+2+2.5} 


3 


e 3 G {0,0.1,0.2,0.3,0.4,0.5,0.6} 


4 


i € {0,^/8,71-/4,371-/8,7172} 


5 


r 3 G {312.5, 625, 1250, 2500, 5000, 10000} 
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FIG. 9. Planar hierarchical system. Comparison between 
the perturbed binary and the unperturbed one. The light 
grey region is the mass quadrupole MQ contribution to the 
waveform rh+, which is not resolved since there are 4070 
orbits. The dark region inside the light grey one is the 
mass octupole plus the current quadrupole MO+CQ contri- 
bution to the waveform rh + . The vertical lines mark the 
time when the simulations are stopped for initial separation 
r 3 G {130, 140, 150, 160, 170} of the inner binary. 



cases the plot shows in light grey the mass quadrupole. 
The waveform looks like a shadow region because com- 
pared to the timescale of the entire evolution a single 
cycle looks like a very high frequency wave. In Fig. |9j 
the binary has completed 3448 orbits, while the inner bi- 
nary of the triple system has completed 4071 orbits and 
the outer binary has completed 7.5 orbits. The mass 
octupole plus the current quadrupole MO+CQ are the 
dark region. Notice that in the triple system MO+CQ 
is modulated by the period of the third body (one cycle 
of modulation corresponds to half an orbit of the third 
body). The perturbation furthermore affects the merger 
time; for the triple system it takes more time for the inner 
binary to merge. We run the simulation for 5 initial inner 
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binary separation r b G {130, 140, 150, 160, 170}. In Fig. [9] 
we mark with vertical lines the time at which the simu- 
lations are stopped. The relative change in the merger 
time 



^3BH — ^2BH 
*2BH 



0.270 ± 0.0025, 



(37) 



is almost constant for this simulations (the standard de- 
viation is 0.0025). We did not observe any particular 
differences in the waveform when changing r b . 



2. Post-Newtonian corrections 

In addition to the comparison to the nonperturbed bi- 
nary system, we use the planar configuration to explore 
the influence of the conservative post-Newtonian correc- 
tions. As in the previous case with initial binary separa- 
tion r b — 150 (which we will denote as full 2.5 PN case), 
we solve the system for equations of motion where we re- 
move the 2 PN part of the Hamiltonian (radiative 1 PN) 
and where we remove both 1 and 2 PN corrections (ra- 
diative Newtonian) . The full 2.5 PN case does not show a 
big difference compared to the radiative 1 PN case. The 
merger phase time changes from t — 4.8372 x 10 7 in the 
first case to t = 4.8132 x 10 7 in the second one. The wave- 
form does not suffer a noticeable change (see Fig.[l0|). On 
the other hand, in the radiative Newtonian case the re- 
sult changes significantly. The merger phase time starts 
later than in previous cases (around t = 5.6388 x 10 7 ). 
For this configuration dynamic which include the radia- 
tive 1 PN corrections seems to be a good approximation. 
However, for the rest of the simulations we employ the 
full 2.5 PN corrections. 



3. Variation of the eccentricity of the external binary 

We analyzed the variation of the waveform as a func- 
tion of the eccentricity of the external binary e%. We 
ran simulations for e3 G {0, 0.1, ... , 0.6}. In this case the 
response to the variation of the eccentricity is better re- 
flected in the combination of i^mli an d hj^£ 3 . Figure 111] 
shows the sum of I = 2, m = 1 and I = m = 3 modes 
(which are the leading components of current quadrupole 
and mass octupole, respectively). The modulation of 
the modes shows two characteristic low- frequencies. If 
we divide the orbit of the external binary in two parts, 
one defined by the true anomaljj^] tp running from tp — 
— 7r/2 to tp = tt/2 and the other by the complement 
tp G [tt/2, 37r/2], it is possible to associate the charac- 
teristic frequencies to each part of the trajectory. We 



1 The true anomaly is defined as the angle which connects the 
periapsis, the main focus and the trajectory of the reduced body. 
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FIG. 10. Successive changes in the waveform due to post- 
Newtonian corrections. Waveform of a radiative Newtonian 
system (bottom), radiative 1 PN system (middle), and full 
2.5 PN system (top). The waveform includes the current and 
mass quadrupole and the mass octupole contributions. The 
vertical dash line at t = 3.81 x 10 7 mark the time when the 
non-perturbed binary system enters the merger phase (see 
Fig.§. 



compute the envelope of the absolute value of the sig- 
nal using a low-pass filter (dark line in Figure 11) for 



2/3 of the total signal (that part of the signal was eas- 
ier to process for high eccentricity). Using the result- 
ing function we compute numerically the local minima. 
The differences between minima are associated with the 
characteristic frequencies. An alternative way to extract 
the characteristic frequencies is by looking at the Fourier 
spectra of the filtered waveform. 

We label the period for tp e [— 7r/2,7r/2] as Ai ap and 
the period for tp G [7r/2,37r/2] as At pcr (at tp = the 
external binary reaches the periapsis and at tp = n the 
apoapsis). Table [TT| shows the results, where we include 
the quotient. 

In the Newtonian case it is possible to compute Ai ap 
and Ai pcr using the conservation of the angular momen- 
tum I and the equation of the orbit (see e.g. |57J. The 
result is 



£3 r3ir/2 

Aiper = — / (1 + ecostp)~ 2 dtp, (38) 

fJ- Jir/2 



;3 r /2 



(1 + e cos tp) 2 dtp, 



(39) 



A* J -tt/2 

where \i is the reduced mass of the binary. The quotient 
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FIG. 11. Sum of modes I — 2, m = 1 and I = m — 3 as func- 
tion of the eccentricity e^. From bottom to top the variation 
of the modes for e e € {0,0.1, . . . ,0.6}. The dark line is the 
envelope of the function for 2/3 of the total simulation. The 
+ marks are the local minima of the envelope. 



TABLE II. Periods At per and At ap and its quotient. The val- 
ues are computed using the averages of the differences between 
the minima (see Figure 11 1 and the errors by the standard de- 
viation. 





0.1 
0.2 
0.3 
0.4 
0.5 
0.6 



Ai pcr [xl0 b 



At ap [xlO 6 ] 



3.1473 
2.3812 
1.7890 
1.3260 
0.9590 
0.6690 
0.4390 



± 0.00020 
± 0.00092 
±0.00180 
± 0.00170 
±0.00160 
±0.00110 
± 0.00330 



3.1427 ± 
3.0700 ± 
2.9950 ± 
2.9160 ± 
2.8370 ± 
2.7530 ± 
2.6670 ± 



0.00053 
0.00120 
0.00051 
0.00110 
0.00110 
0.00170 
0.00400 



0.9985 ± 
1.2890 ± 
1.6750 ± 
2.2000 ± 
2.9580 ± 
4.1180 ± 
6.0700 ± 



0.00023 
0.00100 
0.00190 
0.00360 
0.00610 
0.00920 
0.05500 



between the periods is related to the eccentricity by 
At 



At 



^per " 

a P 2arctan,/iis 

V l+e 



1. 



(40) 
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FIG. 12. Ratio of the periods At pcl and At ap as function of 
the eccentricity. The solid line shows the Newtonian relation. 



sented in Table |TT] and the Newtonian expression ( 40 1 . 
For this case the Newtonian expression represents very 
well the functional behavior of our simulation. 



4- Variation of the inclination angle 

The period of modulation of the I = 3 modes of the 
waveform are related to the period of the third body. On 
the other hand, the amplitude of the I = 3 spherical com- 
ponents of the waveform encode information about the 
inclination angle i. We run simulations with the same 
initial configuration for i E {0, ir/8, 7r/4, 37r/8, 7t/2}. Fig- 
ure 13 shows the variation of the amplitude for the real 
part of the modes h^l 



and h^j^g as a function of i. 



Since the real and the imaginary part of the modes show 
the same behavior, for simplicity we present only the 
analysis of the real part. The real part of h^l 2 is zero 
for planar motion i = 0. However, the contribution of 
this mode increases with i. On the other hand, the con- 
tribution of Re{h ^=3} is maximal in the planar case and 
decreases when i increases. This behavior is symmetric 
with respect to i = n/2 and periodic with period n. 

We estimate the contribution of each mode calculating 
the area which is covered by the real part of the mode, 



A l m (r):=- / |Re{hL(f)}|df, 



Figure 12 shows a comparison between the data pre- 



(41) 



where tf = 4.8372 x 10 7 is the final time of the evolution 
and t — tf — t. We integrate backward in time start- 
ing with the beginning of the merger phase at tf. We 
compute A l m (t) for 8 uniformly spaced times during the 
simulation. We normalize the results using the maximum 
value A max = A l ^ 2 - We denote the normalized area by 
A l m . As an example we show the results for r = in 
Table |III| where we present the relevant modes. In total 
we compute 8 tables similar to the previous one, however 
for brevity we do not present them here. Notice that 
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TABLE III. Variation of A m as a function of the inclination 
angle i. 





— n 


i = 


i = tt/8 


i = 7r/4 


i = 2tt/8 


i = tt/2 


1 


m 


A m 


o 

Z 


u 


0.0019 


0.0019 


0.0021 


0.0024 


0.0026 


2 


1 


0.0000 


0.0007 


0.0013 


0.0016 


0.0018 


2 


2 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


3 





0.0000 


0.0006 


0.0021 


0.0012 


0.0544 


3 


1 


0.0546 


0.0527 


0.0588 


0.0397 


0.1160 


3 


2 


0.0000 


0.0429 


0.0799 


0.1033 


0.1583 


3 


3 


0.2128 


0.2052 


0.1957 


0.1552 


0.2376 
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FIG. 13. Variation of the amplitude of I = 3, m = 2, 3, 
modes as a function of the inclination angle i. Superposition 
of Re{h ^=2} an d Re{h ^=3} as function of i. 
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FIG. 14. Variation of A l m as a function of i for t = tf (upper 
panel) and t = tf/2 (lower panel). 
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the contribution of the I = 2 modes is almost constant 
with respect to the inclination angle i. In Figure 14 we 
show the variation of A^l 2 an d ^m=3 f° r two integra- 
tion times, r = and r = t//2. 

We found that the variation of is wen repre- 

sented by 

Al = 3 

On the other hand, Al =3 



=2 (M) = a(r)|sini|. 



(42) 



\l=3 



(*.*) 



- 3 is well modeled by 
b(r) + c(r) I cos i\ 3/2 , 



(43) 



where the fitting coefficients a, b and c depend on the 
interval of integration. Table IV shows the fitting co- 
efficients as a function of the integration time r. From 
this data it is possible to fit a function to establish the 
functional behavior of the coefficients with respect to the 



integration time. The result is shown in Figure 15 The 
coefficients a, b and c are well represented by 



a(r) 


-T°2 

= one 


(44) 


b(r) 




(45) 


c(r) 


= 7ie 


(46) 



FIG. 15. Functional behavior of the fitting coefficients. The 
coefficients are well described by an exponential decay func- 
tion in t. 



where 





= 8.94 ±0.018, 


(47) 




= (8.352 ± 0.0029) x 10~ 2 , 


(48) 




= 10.17 ±0.21, 


(49) 


& 


= (8.26 ± 0.032) x 10~ 2 , 


(50) 


7i 


= 5.90 ± 0.033, 


(51) 


72 


= (8.305 ± 0.0084) x 10~ 2 . 


(52) 



The asymptotic behavior of the coefficients suggests that 
for long integration times it is possible to consider them 
as constants. 

Alternatively, it is possible to relate the inclination an- 
gle i with the maximum of the modes I = 3, m = 2 and 



I = 2, m = 1. As in Sec. HIB 3 we compute the envelope 
of the modes using a low-pass filter. The upper panel in 
Figure 16 shows the result for the angle i = 7r/4. The 
quotient of the envelope of the modes I = 3, m = 2 and 
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TABLE IV. Fitting coefficients of Eqns. Q and Q. For 
the 8 time intervals we compute the fitting coefficients a, b 
and c. We include the error of each coefficient. 



TABLE V. The maximum of ( 53 1 as a function of the inclina- 



T [Xl0 7 ] 


a(r) [xl0~ 2 ] 


fe(r) [xl0~ 2 ] 


c(t) [x10~ 2 ] 


0.6047 


22.45 ±0.066 


26.99 ±0.116 


15.21 ±0.171 


1.2093 


18.04 ±0.035 


21.57 ±0.062 


12.28 ±0.092 


1.8140 


15.76 ±0.025 


18.90 ± 0.046 


10.75 ± 0.067 


2.4186 


14.28 ±0.019 


17.19 ±0.036 


9.76 ±0.054 


3.0233 


13.21 ±0.016 


15.99 ±0.031 


9.03 ± 0.046 


3.6279 


12.38 ±0.014 


15.07 ±0.027 


8.47 ±0.041 


4.2326 


11.72 ±0.012 


14.33 ±0.025 


8.00 ±0.037 


4.8372 


11.18 ±0.011 


13.73 ±0.024 


7.62 ±0.035 



tion angle i. Listed is the average value of the maxima, while 
the error is given by the standard deviation of the data. 





MaxLRfh l = 3 o h i= Ml 


Variation (%) 


o 


o 





tt/16 


0.1608 ±0.00077 


0.48 


tt/8 


0.335 ±0.0012 


0.36 


3tt/16 


0.538 ± 0.0029 


0.54 


tt/4 


0.806 ± 0.0049 


0.61 


5tt/16 


1.193 ±0.0056 


0.47 


3tt/8 


1.864 ±0.0093 


0.50 


7tt/16 


3.41 ±0.026 


0.75 


tt/2 


5.57 ±0.033 
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FIG. 16. Variation of modes I — 3,m = 2 and I = 2, m = 1 
as a function of the inclination angle i. The upper panel 
shows for i — 7r/4 the absolute value of the modes and its en- 
velope. The lower panel shows the quotient of the envelopes 
I — 3, m = 2 and I = 2, m = 1. Notice that the resulting func- 
tion is almost periodic and does not show the characteristic 
growth close to the merger phase. 



I = 2, m = 1 gives a periodic function which removes the 
growth of the modes close to the merger time. We define 
the function R which rectifies the envelopes as 



-R(fi'n=2i n m=l) 



1=2 



Env[Re{h^ 2 }] 
Env[Re{h^ i} ]' 



(53) 



The lower panel of Figure [l6| shows the result of applying 
(531 to our data. Notice that in the case of i = tt/2 the 
values after t = 3 x 10 7 are a little erratic. For our 
analysis we consider for i = 7r/2 only the points before 
( = 3x 10 7 . 

From the resulting function we compute numerically 
the local maxima of d53|. Table |V] shows the result. 



For this purpose we perform additional simulations for 
angles 7r/16, 37r/16, 57r/16 and 77r/16. We fit to the data 
the function f{i) = aie b l , where a = 0.65 ± 0.034 and 
b = 0.69 ± 0.024. Figure 17 shows the result, notice that 



4.8 



3.6 



1.2 



h 1 1 1 1 1 1 r 

• — i — ■ data , 

aie bi 

a = 0.65 ± 0.034 
b = 0.69 ± 0.024 



A 



1~ 



_i i i i i i i i_ 



tt/8 n/4 3tt/8 tt/2 

i 



FIG. 17. The maximum of (531 as function of the inclination 



angle i. The functional behavior is well represented by the 
function aie * . 



the functional behavior is well represented by the fitted 
function. 

In both cases, using the relative "area" of the modes or 
the maximum of the "rectified" modes, we obtain quite 
a simple behavior. The advantage of the second method 
is that it does not depend on the integration time r. 



5. Initial separation of the external binary 

The last numerical experiment examines the depen- 
dence on the initial separation of the external binary r%. 
We set the value of r 3 to 312.5, 625, 1250, 2500, 5000, and 
10000. For r 3 = 312.5 the external body is ejected from 
the binary after a few orbits, the other configurations are 
stable. 



Figure 18 shows the sum of the mass octupole and 
current quadrupolc contributions to the waveform. The 
frequency of the modulation of the waveform increases 
when the separation and hence the orbital period of the 
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FIG. 18. Planar hierarchical system. Modulation of the mass 
octupole plus the current quadrupole as function of the initial 
separation of the external body. The initial separation rz 
takes the values 625 (a), 1250 (b), 2500 (c), 5000 (d) and 
10000 (e). Shown on the left is the evolution for t G [0, 10 7 ] 
and on the right for t G [4 x 10 7 , 5 x 10 7 ]. 



external binary is decreased. One orbit of the external bi- 
nary corresponds to the time between two of the nodes of 
the mass octupole plus current quadrupole contribution 



shown in Figure 18 The influence of a third body is not 
clearly defined when the period of the external binary is 
similar to the inner binary. For small separations, on the 
scale shown there is no modulation of the waves visible 
(see Figure [l8| (a) and (b)). When the initial separation 
of the external binary is increased, at some distance most 
of the inspiral and merger of the inner binary happens 
before the external binary completes one orbit. 



IV. DISCUSSION 

We performed post-Newtonian simulations for a se- 
lection of hierarchical configurations as an example for 
a three-body system, and we analyzed the waveforms. 
Based on these simulations we examined a number of 
different physical aspects of the system. 

First of all, looking at the mass octupole and current 
quadrupole part of the waveform, it is possible to distin- 
guish between such a hierarchical (also called Jacobian) 
triple system and a binary system, an issue that has been 
discussed in [2J |3] . 

In terms of the merger time, the perturbed binary 
merges later. For mass ratio 10:20:1, the delay of the 
merger is 27% compared to the binary with 10:20, which 
is perhaps surprisingly large. However, let us note that 



even a small perturbation due to a third object can have 
a large effect when integrated over about 4000 orbits of 
the inner binary (i.e. there is less than a 0.01% delay 
per orbit). As we have shown, the delay depends only 
very weakly on the inclination angle or the distance to 
the third body, see Figs. [13] and [l8j This may be ex- 
pected since the force due to the third body periodically 
increases but also decreases the force between the objects 
of the inner binary (depending on the orientation of the 
binary with respect to the third body) , which apparently 
averages out over several orbits of the inner binary. As 
a cross check we also performed simulations where the 
third mass approaches zero, and in this case the merger 
time does approach that of the binary. 

As far as the approximation method is concerned, we 
find that there is a significant difference in the merger 
time for a system which includes Newtonian dynamics 
and 2.5 PN radiation compared to the inclusion of 1 PN 
or 2 PN corrections to the dynamics. The inclusion of 
1 PN corrections to the conservative part of the Hamilto- 
nian produces a change of 16% in the merger time. How- 
ever, the inclusion of 2 PN corrections does not make a 
significant difference to either the waveform or the merger 
time (only around 0.5%). 

The variation of the eccentricity of the external binary 
shows that the period of the third body is well described 
by the Newtonian dynamics. From the modulation of 
the waveform modes (particularly from the sum of the 
I = 2, to = 1 and I = m = 3 modes), it is possible 
to distinguish two frequencies which are related to the 
eccentricity via a Newtonian expression. 

We established a link between the amplitude of the 
I = 3, m = 2 and I = to = 3 modes and the angle of the 
osculating orbital planes. In order to extract the infor- 
mation given by the waves we used two methods. First, 
we used the relative area covered by the I — 3, m = 2 
and / = to = 3 modes with respect to the area covered 
by the mode I = m = 2. In this case the contribution 
of the I = 3, to = 2 mode is particularly simple. It is 
zero for planar motion and increases as a sine function 
of the inclination angle. The second method is based on 
the quotient of the envelope of the I — 3, to = 2 mode 
and the envelope of the I = 2, m = 1 mode. The re- 
sulting function is almost periodic and does not contain 
the characteristic growth of the waveforms close to the 
merger phase. In this case, it is possible to relate the 
inclination angle to the amplitude of the resulting func- 
tion. The modulation produced by the third body on 
the I — 3 modes characterizes the period of the external 
binary. Decreasing the initial separation of the external 
body produces a higher frequency modulation, until it is 
no longer possible to discern a well defined modulation 
of the waveform. In our simulations, when there are no 
well defined internal and external binaries the system is 
not stable. 

Our results provide additional evidence to a conjecture 
first stated in [5] , that in order to characterize a system of 
n compact objects, it is necessary to perform an analysis 
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of the waveform which includes at least the I < n modes. 
As we showed in the last numerical experiment, when the 
third body is close to the binary it is not evident how to 
extract information related to the dynamics of a partic- 
ular body. It is necessary to perform a detailed study of 
nonhierarchical triple systems to determine how much in- 
formation we can extract from more general cases. More 
detailed statements based on the higher modes of the 
waveform are possible but require an extensive param- 
eter study. Other configurations include for example a 
massive compact object perturbing a binary, or the scat- 
tering and capture of a third body. The present examples 
showed the type of characterization that are possible with 
the techniques developed above. 



As a final comment, let us point out that chaotic be- 
havior of triple systems is well known in the Newtonian 
case (see e.g. [I] and references therein). For binaries, 
it is known that chaos appears when using certain post- 
Newtonian approximations for systems of spinning bina- 
ries (see e.g. As a natural generalization of the 
Newtonian case we expect that the three-body problem 
exhibits chaotic behavior as well. An important ques- 
tion is, how does the emission of gravitational radiation 
change the chaotic properties of the system? We consider 
this a topic for future study. 
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Appendix A: First and second post-Newtonian 
Hamiltonian 



Here we reproduce in our notation the Hamiltonian 
given in [34] , with some factorizations and changes in the 
summation of the terms Tl and T2, which are marked by 
braces below. Our version (worked out with G. Schafer) 
fixes the typos noted in [35] , giving a formula equivalent 
to [35j but written in a different way. The issue is how the 
four-point functions of [37] are reduced to explicit triple 
sums for a three-body problem. The first and second 
post-Newtonian Hamiltonians are 



4 r ab \ m a 

a b^a 



a \ a/ a b _ ia 

7p a ■ Pb ~ (fl a b ■ Pa)(n a b ' Pb)j 
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Appendix B: Lagrange triangle solution waveform 

Here we summarize the expressions for the mass 
quadrupole, mass octupole, and current quadrupole 
waveforms for each polarization of the Lagrange trian- 
gle solution. See 3J for details on the calculation of this 
expression. We denote by a := r 12 = r 13 = 7' 2 3 the sepa- 
ration between each pair of bodies, toi, TO2 and TO3 are 



the dimensionless mass parameters, uj = a~ 3 / 2 is the or- 
bital frequency, r is the distance from the observer to the 
source and 9 is the observational direction. We define the 
following auxiliary quantities: 



Mi : = + TOj-TOfc + TO 2 , (Bl) 

0i := 0, (B2) 
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t>2 '■= arccos 



arccos 



Mi + M3 - 1 
2/xi/i 3 

Mi + mI - 1 

2/xiM2 



(B3) 
(B4) 



where j ^ i, k ^ The plus and cross polarizations 
of the mass quadrupole waveform are 



rh!f Q = -(3 + cos 26)a 2 uj 2 V m^- cos(2(wi + (pi)) 
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rh 1 ^® = —4 cos 9a 2 uj 2 nii/j- sin(2(wi 
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fa)), (B6) 



the expressions for the current quadrupole are 
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sin 9 TOj/zf cos(w< + < 



(B7) 



eg 



2a 3 w 3 



sin(20) J^m^ 3 sin(wi + <fo ), (B8) 



»=i 



and the waveforms for the mass octupole are given by 
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